We present recent work on efficient numerical methods for inverse problems governed by transport equations. The control variable of our problem is either a stationary or a time-dependent velocity field. The variational problem formulation includes of a data mismatch term, a Sobolev norm to regularize the control, and mixed-type partial differential equations (PDEs) as constraints. In its simplest form, the constraints consist of a hyperbolic transport equation. We extend this formulation by introducing additional constraints on the divergence of the velocity field to achieve incompressible or near-incompressible flows. Alternatively, we introduce a geodesic equation on the group of diffeomorphisms, resulting in an initial value control problem. In this case, the control variable becomes an initial momentum or initial velocity field at time t = 0. The underlying inverse problem is highly non-linear, inherently ill-posed, and infinite-dimensional in the continuum, resulting in high-dimensional, ill-conditioned inversion operators after discretization. This represents significant mathematical and computational challenges. We employ a Gauss-Newton-Krylov method for numerical optimization, utilizing spectral discretization in space and a semi-Lagrangian method for numerical time integration. Our code has been deployed in modern heterogeneous high-performance computing platforms. We have designed dedicated mixed-precision computational kernels to maximize throughput and efficiency. We discuss effective numerical methods for evaluating forward and adjoint operators, preconditioning the reduced space Hessian, effective strategies for numerical optimization, and the scalability of our solver. We test the performance of the proposed numerical scheme on synthetic and real-world data.